Bayesian GLM Part7

Author

Murray Logan

Published

08/07/2025

1 Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(HDInterval) # for HPD intervals
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(posterior) # for posterior draws
library(ggeffects) # for partial effects plots
library(patchwork) # for multi-panel figures
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # framework for stats, modelling and visualisation
library(geoR) # framework for stats, modelling and visualisation
library(modelsummary) # for data and model summaries
theme_set(theme_grey()) # put the default ggplot theme back
source("helperFunctions.R")

2 Scenario

Mitchell et al. (2019) presented a data set in which male guppy fish where placed individually into long, narrow and shallow tanks that predominately only permitted the fish to swim in two dimensions. The activity rate (cumulative distance swam) of each guppy fish was remotely recorded via tracking software during 20 minute trials conducted each day for 14 days. The order in which individual guppy fish were recorded each day was randomised and the time of day of each trial was recorded.

The researchers were interested in whether the the guppy fish become habituated to the experimental conditions and altered their activity patterns over time.

The data are in the file mitchell.csv in the data folder.

Table 1: Format of the mitchell.csv data file

ID	TOD	Day	Dist_moved
A39	0.39	7	3439.93
A3	0.39	7	2795.08
A13	0.39	7	5258.79
...	...

Table 2: Description of the variables in the mitchell data file

ID:	Individual guppy fish ID - Random variable
TOD:	Time of day of observation as a proportion - Predictor variable
Day:	Day post entry into tank - Predictor variable
Dist_moved:	Total distance (cm) moved in 20 min trial - Response variable

3 Read in the data

mitchell <- read_csv("../data/mitchell.csv", trim_ws = TRUE)

Rows: 1626 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): ID, Date, StartDate
dbl  (5): Mass, Day, TOD, Dist_moved, Dist
time (1): TimeOfDay

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(mitchell)

Rows: 1,626
Columns: 9
$ ID         <chr> "A39", "A3", "A13", "A23", "A59", "A69", "A41", "A9", "A35"…
$ TimeOfDay  <time> 09:28:32, 09:28:32, 09:28:32, 09:28:32, 09:28:32, 09:28:32…
$ Date       <chr> "6/03/2018", "6/03/2018", "6/03/2018", "6/03/2018", "6/03/2…
$ StartDate  <chr> "27-Feb-18", "27-Feb-18", "27-Feb-18", "27-Feb-18", "27-Feb…
$ Mass       <dbl> 0.066, 0.072, 0.105, 0.080, 0.077, 0.092, 0.073, 0.075, 0.0…
$ Day        <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,…
$ TOD        <dbl> 0.39, 0.39, 0.39, 0.39, 0.39, 0.39, 0.39, 0.42, 0.42, 0.42,…
$ Dist_moved <dbl> 3439.930, 2795.080, 5258.790, 3402.740, 7672.660, 2400.660,…
$ Dist       <dbl> 58.65092, 52.86852, 72.51752, 58.33301, 87.59372, 48.99653,…

## Explore the first 6 rows of the data
head(mitchell)

# A tibble: 6 × 9
  ID    TimeOfDay Date      StartDate  Mass   Day   TOD Dist_moved  Dist
  <chr> <time>    <chr>     <chr>     <dbl> <dbl> <dbl>      <dbl> <dbl>
1 A39   09:28:32  6/03/2018 27-Feb-18 0.066     7  0.39      3440.  58.7
2 A3    09:28:32  6/03/2018 27-Feb-18 0.072     7  0.39      2795.  52.9
3 A13   09:28:32  6/03/2018 27-Feb-18 0.105     7  0.39      5259.  72.5
4 A23   09:28:32  6/03/2018 27-Feb-18 0.08      7  0.39      3403.  58.3
5 A59   09:28:32  6/03/2018 27-Feb-18 0.077     7  0.39      7673.  87.6
6 A69   09:28:32  6/03/2018 27-Feb-18 0.092     7  0.39      2401.  49.0

str(mitchell)

spc_tbl_ [1,626 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ ID        : chr [1:1626] "A39" "A3" "A13" "A23" ...
 $ TimeOfDay : 'hms' num [1:1626] 09:28:32 09:28:32 09:28:32 09:28:32 ...
  ..- attr(*, "units")= chr "secs"
 $ Date      : chr [1:1626] "6/03/2018" "6/03/2018" "6/03/2018" "6/03/2018" ...
 $ StartDate : chr [1:1626] "27-Feb-18" "27-Feb-18" "27-Feb-18" "27-Feb-18" ...
 $ Mass      : num [1:1626] 0.066 0.072 0.105 0.08 0.077 0.092 0.073 0.075 0.09 0.077 ...
 $ Day       : num [1:1626] 7 7 7 7 7 7 7 7 7 7 ...
 $ TOD       : num [1:1626] 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.42 0.42 0.42 ...
 $ Dist_moved: num [1:1626] 3440 2795 5259 3403 7673 ...
 $ Dist      : num [1:1626] 58.7 52.9 72.5 58.3 87.6 ...
 - attr(*, "spec")=
  .. cols(
  ..   ID = col_character(),
  ..   TimeOfDay = col_time(format = ""),
  ..   Date = col_character(),
  ..   StartDate = col_character(),
  ..   Mass = col_double(),
  ..   Day = col_double(),
  ..   TOD = col_double(),
  ..   Dist_moved = col_double(),
  ..   Dist = col_double()
  .. )
 - attr(*, "problems")=<externalptr>

mitchell |> datawizard::data_codebook()

mitchell (1626 rows and 9 variables, 9 shown)

ID | Name       | Type      | Missings |            Values |           N
---+------------+-----------+----------+-------------------+------------
1  | ID         | character | 0 (0.0%) |                A1 |  14 ( 0.9%)
   |            |           |          |               A11 |  14 ( 0.9%)
   |            |           |          |               A13 |  14 ( 0.9%)
   |            |           |          |               A15 |  14 ( 0.9%)
   |            |           |          |               A17 |  14 ( 0.9%)
   |            |           |          |               A19 |  14 ( 0.9%)
   |            |           |          |               A21 |  14 ( 0.9%)
   |            |           |          |               A23 |  14 ( 0.9%)
   |            |           |          |               A25 |  14 ( 0.9%)
   |            |           |          |               A27 |  14 ( 0.9%)
   |            |           |          |             (...) |            
---+------------+-----------+----------+-------------------+------------
2  | TimeOfDay  | numeric   | 0 (0.0%) |          07:57:31 |   7 ( 0.4%)
   |            |           |          |          08:08:00 |   7 ( 0.4%)
   |            |           |          |          08:16:18 |   2 ( 0.1%)
   |            |           |          |          08:17:16 |   8 ( 0.5%)
   |            |           |          |          08:17:56 |   5 ( 0.3%)
   |            |           |          |          08:18:36 |   7 ( 0.4%)
   |            |           |          |          08:20:22 |   6 ( 0.4%)
   |            |           |          |          08:23:28 |   8 ( 0.5%)
   |            |           |          |          08:26:07 |   5 ( 0.3%)
   |            |           |          |          08:30:43 |  10 ( 0.6%)
   |            |           |          |             (...) |            
---+------------+-----------+----------+-------------------+------------
3  | Date       | character | 0 (0.0%) |         1/04/2018 |  40 ( 2.5%)
   |            |           |          |         1/05/2018 |  40 ( 2.5%)
   |            |           |          |        10/03/2018 |  39 ( 2.4%)
   |            |           |          |        10/04/2018 |  40 ( 2.5%)
   |            |           |          |        11/03/2018 |  39 ( 2.4%)
   |            |           |          |        11/04/2018 |  40 ( 2.5%)
   |            |           |          |        12/03/2018 |  39 ( 2.4%)
   |            |           |          |        12/04/2018 |  40 ( 2.5%)
   |            |           |          |        13/03/2018 |  39 ( 2.4%)
   |            |           |          |        14/03/2018 |  39 ( 2.4%)
   |            |           |          |             (...) |            
---+------------+-----------+----------+-------------------+------------
4  | StartDate  | character | 0 (0.0%) |         17-Apr-18 | 560 (34.4%)
   |            |           |          |         23-Mar-18 | 520 (32.0%)
   |            |           |          |         27-Feb-18 | 546 (33.6%)
---+------------+-----------+----------+-------------------+------------
5  | Mass       | numeric   | 0 (0.0%) |      [0.05, 0.12] |        1626
---+------------+-----------+----------+-------------------+------------
6  | Day        | numeric   | 0 (0.0%) |           [7, 20] |        1626
---+------------+-----------+----------+-------------------+------------
7  | TOD        | numeric   | 0 (0.0%) |      [0.37, 0.56] |        1626
---+------------+-----------+----------+-------------------+------------
8  | Dist_moved | numeric   | 0 (0.0%) | [101.13, 11833.9] |        1626
---+------------+-----------+----------+-------------------+------------
9  | Dist       | numeric   | 0 (0.0%) |   [10.06, 108.78] |        1626
------------------------------------------------------------------------

mitchell |> modelsummary::datasummary_skim()

Warning: These variables were omitted because they include more than 50 levels:
ID.

	Unique	Missing Pct.	Mean	SD	Min	Median	Max
Mass	51	0	0.1	0.0	0.0	0.1	0.1
Day	14	0	13.5	4.1	7.0	13.0	20.0
TOD	97	0	0.5	0.0	0.4	0.4	0.6
Dist_moved	1625	0	3917.0	2141.8	101.1	3811.0	11833.9
Dist	1625	0	59.9	18.1	10.1	61.7	108.8
		N	%
Date	1/04/2018	40	2.5
	1/05/2018	40	2.5
	10/03/2018	39	2.4
	10/04/2018	40	2.5
	11/03/2018	39	2.4
	11/04/2018	40	2.5
	12/03/2018	39	2.4
	12/04/2018	40	2.5
	13/03/2018	39	2.4
	14/03/2018	39	2.4
	15/03/2018	39	2.4
	16/03/2018	39	2.4
	17/03/2018	39	2.4
	18/03/2018	39	2.4
	19/03/2018	39	2.4
	2/04/2018	40	2.5
	2/05/2018	40	2.5
	24/04/2018	40	2.5
	25/04/2018	40	2.5
	26/04/2018	40	2.5
	27/04/2018	40	2.5
	28/04/2018	40	2.5
	29/04/2018	40	2.5
	3/04/2018	40	2.5
	3/05/2018	40	2.5
	30/03/2018	40	2.5
	30/04/2018	40	2.5
	31/03/2018	40	2.5
	4/04/2018	40	2.5
	4/05/2018	40	2.5
	5/04/2018	40	2.5
	5/05/2018	40	2.5
	6/03/2018	39	2.4
	6/04/2018	40	2.5
	6/05/2018	40	2.5
	7/03/2018	39	2.4
	7/05/2018	40	2.5
	8/03/2018	39	2.4
	8/04/2018	40	2.5
	9/03/2018	39	2.4
	9/04/2018	40	2.5
StartDate	17-Apr-18	560	34.4
	23-Mar-18	520	32.0
	27-Feb-18	546	33.6

mitchell |> modelsummary::datasummary_skim(by = "Day")

Warning: These variables were omitted because they include more than 50 levels:
ID.

	Day	Unique	Missing Pct.	Mean	SD	Min	Median	Max
Mass	7.0	51	0	0.1	0.0	0.0	0.1	0.1
	8.0	51	0	0.1	0.0	0.0	0.1	0.1
	9.0	51	0	0.1	0.0	0.0	0.1	0.1
	10.0	51	0	0.1	0.0	0.0	0.1	0.1
	11.0	51	0	0.1	0.0	0.0	0.1	0.1
	12.0	51	0	0.1	0.0	0.0	0.1	0.1
	13.0	51	0	0.1	0.0	0.0	0.1	0.1
	14.0	51	0	0.1	0.0	0.0	0.1	0.1
	15.0	45	0	0.1	0.0	0.0	0.1	0.1
	16.0	51	0	0.1	0.0	0.0	0.1	0.1
	17.0	51	0	0.1	0.0	0.0	0.1	0.1
	18.0	51	0	0.1	0.0	0.0	0.1	0.1
	19.0	51	0	0.1	0.0	0.0	0.1	0.1
	20.0	51	0	0.1	0.0	0.0	0.1	0.1
TOD	7.0	12	0	0.5	0.0	0.4	0.5	0.5
	8.0	11	0	0.5	0.1	0.4	0.4	0.6
	9.0	18	0	0.5	0.1	0.4	0.4	0.6
	10.0	18	0	0.5	0.0	0.4	0.5	0.5
	11.0	19	0	0.4	0.0	0.4	0.4	0.5
	12.0	19	0	0.4	0.0	0.4	0.4	0.5
	13.0	20	0	0.5	0.1	0.4	0.4	0.5
	14.0	21	0	0.5	0.1	0.4	0.4	0.5
	15.0	13	0	0.5	0.1	0.4	0.5	0.6
	16.0	17	0	0.5	0.0	0.4	0.4	0.5
	17.0	17	0	0.5	0.0	0.4	0.5	0.6
	18.0	20	0	0.5	0.1	0.4	0.5	0.5
	19.0	17	0	0.5	0.0	0.4	0.4	0.6
	20.0	18	0	0.5	0.0	0.4	0.5	0.5
Dist_moved	7.0	119	0	3766.9	2207.2	101.1	3402.7	9037.0
	8.0	119	0	3609.6	2066.8	269.8	3589.3	8343.5
	9.0	119	0	3468.2	1996.3	251.3	3396.0	8581.9
	10.0	119	0	3525.7	1908.6	472.4	3256.8	8343.1
	11.0	119	0	3603.1	2032.4	341.6	3389.6	8508.7
	12.0	119	0	3682.0	2081.9	577.9	3425.5	9167.5
	13.0	119	0	3918.3	2181.0	415.2	3691.7	8840.6
	14.0	119	0	4018.1	2107.9	566.2	3928.4	8881.7
	15.0	79	0	4031.3	2212.8	685.7	3851.9	9619.0
	16.0	119	0	4199.5	2241.9	438.1	4215.6	10238.8
	17.0	119	0	4234.3	2222.7	549.3	4107.0	9439.1
	18.0	119	0	4315.3	2128.9	563.4	4684.7	9420.1
	19.0	119	0	4238.4	2202.5	656.7	4053.9	9186.6
	20.0	119	0	4266.2	2220.8	433.4	4454.1	11833.9
Dist	7.0	119	0	58.3	19.4	10.1	58.3	95.1
	8.0	119	0	57.3	18.1	16.4	59.9	91.3
	9.0	119	0	56.2	17.7	15.9	58.3	92.6
	10.0	119	0	56.9	16.9	21.7	57.1	91.3
	11.0	119	0	57.3	18.0	18.5	58.2	92.2
	12.0	119	0	58.1	17.6	24.0	58.5	95.7
	13.0	119	0	59.8	18.6	20.4	60.8	94.0
	14.0	119	0	60.9	17.7	23.8	62.7	94.2
	15.0	79	0	60.8	18.4	26.2	62.1	98.1
	16.0	119	0	62.2	18.4	20.9	64.9	101.2
	17.0	119	0	62.6	18.0	23.4	64.1	97.2
	18.0	119	0	63.3	17.6	23.7	68.4	97.1
	19.0	119	0	62.6	18.0	25.6	63.7	95.8
	20.0	119	0	62.8	18.1	20.8	66.7	108.8
		N	%
Date	1/04/2018	40	2.5
	1/05/2018	40	2.5
	10/03/2018	39	2.4
	10/04/2018	40	2.5
	11/03/2018	39	2.4
	11/04/2018	40	2.5
	12/03/2018	39	2.4
	12/04/2018	40	2.5
	13/03/2018	39	2.4
	14/03/2018	39	2.4
	15/03/2018	39	2.4
	16/03/2018	39	2.4
	17/03/2018	39	2.4
	18/03/2018	39	2.4
	19/03/2018	39	2.4
	2/04/2018	40	2.5
	2/05/2018	40	2.5
	24/04/2018	40	2.5
	25/04/2018	40	2.5
	26/04/2018	40	2.5
	27/04/2018	40	2.5
	28/04/2018	40	2.5
	29/04/2018	40	2.5
	3/04/2018	40	2.5
	3/05/2018	40	2.5
	30/03/2018	40	2.5
	30/04/2018	40	2.5
	31/03/2018	40	2.5
	4/04/2018	40	2.5
	4/05/2018	40	2.5
	5/04/2018	40	2.5
	5/05/2018	40	2.5
	6/03/2018	39	2.4
	6/04/2018	40	2.5
	6/05/2018	40	2.5
	7/03/2018	39	2.4
	7/05/2018	40	2.5
	8/03/2018	39	2.4
	8/04/2018	40	2.5
	9/03/2018	39	2.4
	9/04/2018	40	2.5
StartDate	17-Apr-18	560	34.4
	23-Mar-18	520	32.0
	27-Feb-18	546	33.6

4 Data preparation

Let start by declaring the categorical variables and random effect as factors.

mitchell <- mitchell |>
  mutate(
    ID = factor(ID),
    fDay = factor(Day)
  )

5 Exploratory data analysis

As the response represents a strictly positive real number (total distance moved), it would make sense to start by considering a Gamma model.

Note, also that the response variable are relatively large distances (as the units of distance used were centimeters). It might be worth dividing these values by 100 or even 1000 during the modelling process.

Model formula: \[ y_i/1000 \sim{} \mathcal{Gamma}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where $\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$ are vectors of the fixed and random effects parameters respectively and $\bf{X}$ is the model matrix representing the overall intercept and effects of day and time of day on the total distance moved during the trial. $\bf{Z}$ represents a cell means model matrix for the random intercepts associated with individual guppies.

mitchell |> ggplot(aes(y = Dist_moved)) +
  geom_boxplot()

## mitchell |> ggplot(aes(y = Dist)) +
##   geom_boxplot()
mitchell |> ggplot(aes(y = Dist_moved)) +
  geom_boxplot() +
  scale_y_log10()

mitchell |> ggplot(aes(y = TOD)) +
  geom_boxplot() +
  scale_y_log10()

mitchell |> ggplot(aes(y = Day)) +
  geom_boxplot()

## mitchell |> ggplot(aes(y = Dist, x = Day)) +
##   geom_point() +
##   geom_line(aes(group = ID))
mitchell |> ggplot(aes(y = Dist_moved, x = TOD)) +
  geom_point() +
  geom_line(aes(group = ID))

mitchell |> ggplot(aes(y = Dist, x = Day)) +
  geom_point() +
  scale_y_log10() +
  scale_x_log10() +
  geom_smooth(method = "lm") +
  facet_wrap(~ID)

`geom_smooth()` using formula = 'y ~ x'

mitchell |> ggplot(aes(y = Dist, x = TOD)) +
  geom_point() +
  scale_y_log10() +
  scale_x_log10() +
  geom_smooth(method = "lm") +
  facet_wrap(~ID)

`geom_smooth()` using formula = 'y ~ x'

mitchell |> ggplot(aes(y = Day, x = TOD)) +
  geom_point()

6 Fit the model

mitchell_sub <- mitchell |>
  ## filter(ID == "A63") |>
  filter(ID == "C69") |>
  droplevels()

brms

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

mitchell.form <- bf(I(Dist_moved/1000) ~ scale(Day) + scale(TOD),
                family=Gamma(link='log'))
options(width=150)
mitchell.form |> get_prior(data = mitchell_sub)

                  prior     class     coef group resp dpar nlpar lb ub       source
                 (flat)         b                                           default
                 (flat)         b scaleDay                             (vectorized)
                 (flat)         b scaleTOD                             (vectorized)
 student_t(3, 1.3, 2.5) Intercept                                           default
      gamma(0.01, 0.01)     shape                                 0         default

options(width=80)

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

$\beta_0$: normal centred at 134 with a standard deviation of 65
- mean of 134: since median(log(mitchell$Dist_moved/1000))
- sd pf 0.7: since mad(log(mitchell$Dist_moved/1000))
$\beta_1$: normal centred at 0 with a standard deviation of 2
- sd of 1: since mad(log(Dist_moved/1000))/mad(scale(Day))
$\sigma$: half-t centred at 0 with a standard deviation of 65 OR
- sd pf 65: since mad(log(mitchell$Dist_moved/1000))
$\sigma$: gamma with shape parameters of 2 and 1
$\beta_0$: normal centred at 1.5 with a standard deviation of 1.5
- mean of 1.5: since mean(log(mitchell$Dist_moved/1000)) or mean(asinh(mitchell$Dist_moved/(2*1))/log(exp(1)))
- sd of 1.5: since sd(log(mitchell$Dist_moved/1000)) or sd(asinh(mitchell$Dist_moved/(2*1))/log(exp(1)))

mitchell_sub |>
  mutate(Dist = Dist) |>
  summarise(
    Int_mu = median(Dist),
    Int_sd = mad(Dist),
    b1_sd = mad(Dist) / mad(scale(Day)),
    b2_sd = mad(Dist) / mad(scale(TOD))
  )

# A tibble: 1 × 4
  Int_mu Int_sd b1_sd b2_sd
   <dbl>  <dbl> <dbl> <dbl>
1   61.5   13.4  10.8  9.97

## mitchell |>
##   mutate(lDist = log(Dist_moved/1)) |>
##   summarise(
##     Int_mu = median(lDist),
##     Int_sd = mad(lDist),
##     b1_sd = mad(lDist)/mad(scale(Day)),
##     b2_sd = mad(lDist)/mad(scale(TOD))
##     )
## standist::visualize("normal(1.4, 1)", xlim=c(0,20))
## standist::visualize("student_t(3, 0, 1)",
##                     xlim=c(-10,25))

## priors <- prior(normal(1.4, 1), class = 'Intercept') +
##     prior(normal(0, 1), class = 'b') +
##     prior(student_t(3, 0, 1), class = 'sd')
##     prior(gamma(2, 1), class = 'shape')

priors <- prior(normal(60, 25), class = "Intercept") +
  prior(normal(0, 15), class = "b") +
  ## prior(student_t(3, 0, 25), class = 'sd') +
  prior(student_t(3, 0, 25), class = "sigma")
## prior(gamma(2, 1), class = 'shape')
## mitchell.form <- bf(Dist_moved ~ scale(Day) + scale(TOD) + (scale(Day) + scale(TOD)|ID),
##                 ## family=Gamma(link='log'))
##                 family=gaussian())
mitchell.form <- bf(Dist ~ scale(Day) + scale(TOD),
  ## family=Gamma(link='log'))
  family = gaussian()
)
mitchell.brm2 <- brm(mitchell.form,
  data = mitchell_sub,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 100,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

mitchell.brm2 |>
  conditional_effects("Day") |>
  plot(points = TRUE)

mitchell.brm2 |>
  conditional_effects("Day") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm2 |>
  conditional_effects("TOD") |>
  plot(points = TRUE)

mitchell.brm2 |>
  conditional_effects("TOD") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

The above seem sufficiently wide whilst at the same time not providing any encouragement for the sampler to wander off into very unsupported areas.

mitchell.brm3 <- update(mitchell.brm2,
  sample_prior = "yes",
  chains = 3, cores = 3,
  refresh = 100
)

The desired updates require recompiling the model

Compiling Stan program...

Start sampling

mitchell.brm3 |>
  conditional_effects("Day") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm3 |>
  conditional_effects("TOD") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm3 |> SUYR_prior_and_posterior()

7 MCMC sampling diagnostics

brms

stan plots

The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

mitchell.brm3$fit |> stan_trace()

mitchell.brm3$fit |> stan_trace(inc_warmup = TRUE)

The chains appear well mixed and very similar

stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

mitchell.brm3$fit |> stan_ac()

There is no evidence of auto-correlation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

mitchell.brm3$fit |> stan_rhat()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

mitchell.brm3$fit |> stan_ess()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

mitchell.brm3$fit |> stan_dens(separate_chains = TRUE)

8 Model validation

brms

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

mitchell.brm3 |> pp_check(type = "dens_overlay", ndraws = 100)

The model draws appear to represent the shape of the observed data reasonably well

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

mitchell.brm3 |> pp_check(type = "error_scatter_avg")

Using all posterior draws for ppc type 'error_scatter_avg' by default.

The predictive error seems to be related to the predictor - the model performs poorest at higher recruitments.

intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

mitchell.brm3 |> pp_check(type = "intervals")

Using all posterior draws for ppc type 'intervals' by default.

## mitchell.brm3 |> pp_check(group='DENSITY', type='intervals')

Whilst the modelled predictions do a reasonable job of representing the observed data, the observed data do appear to be more varied than the model is representing.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(mitchell.brm3)

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

preds <- mitchell.brm3 |> posterior_predict(nsamples = 250, summary = FALSE)

Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.

## mitchell.resids <- createDHARMa(simulatedResponse = t(preds),
##                             observedResponse = mitchell$Dist_moved/1000,
##                             fittedPredictedResponse = apply(preds, 2, median),
##                             integerResponse = TRUE)
mitchell.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = mitchell_sub$Dist,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
mitchell.resids |> plot()

mitchell.resids |> testDispersion()


    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.49418, p-value = 0.416
alternative hypothesis: two.sided

mitchell.resids <- make_brms_dharma_res(mitchell.brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(mitchell.resids)) +
  wrap_elements(~ plotResiduals(mitchell.resids, form = factor(rep(1, nrow(mitchell_sub))))) +
  wrap_elements(~ plotResiduals(mitchell.resids, quantreg = TRUE)) +
  wrap_elements(~ testDispersion(mitchell.resids))

mitchell.resids |> testTemporalAutocorrelation(time = mitchell_sub$Day)


    Durbin-Watson test

data:  simulationOutput$scaledResiduals ~ 1
DW = 2.0389, p-value = 0.9405
alternative hypothesis: true autocorrelation is not 0

## mitchell.resid1 <- mitchell.resids |>
##   recalculateResiduals(group = mitchell$Day, aggregateBy = mean)
##   ## recalculateResiduals(group = interaction(mitchell$Day,  mitchell$ID),  aggregateBy = mean)
## mitchell.resid1 |> testTemporalAutocorrelation(time=unique(mitchell$Day))
autocor_check(mitchell_sub, mitchell.brm3, variable = "Day", n.sim = 250)

variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram

Using  as id variables

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

`geom_smooth()` using formula = 'y ~ x'

residuals(mitchell.brm3)[, "Estimate"] |> acf()

Conclusions:

the simulated residuals do suggest that there might be a dispersion issue
it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.

9 Refit model

## prior(gamma(2, 1), class = 'shape')
priors <- prior(normal(60, 10), class = "Intercept") +
  prior(normal(0, 10), class = "b") +
  ## prior(student_t(3, 0, 10), class = 'sd')  +
  prior(student_t(3, 0, 10), class = "sigma") #+
## prior(lkj_corr_cholesky(1), class = 'cor')
## prior(gamma(0.01, 0.01), class = 'shape')

mitchell.form <- bf(
  Dist ~ scale(Day) + scale(TOD) +
    ar(time = Day, gr = ID, p = 1, cov = TRUE),
  family = gaussian()
)
## family=lognormal())
## family=Gamma(link='log'))
get_prior(mitchell.form, data = mitchell_sub)

                    prior     class     coef group resp dpar nlpar lb ub
                   (flat)        ar                                -1  1
                   (flat)         b                                     
                   (flat)         b scaleDay                            
                   (flat)         b scaleTOD                            
 student_t(3, 61.5, 13.4) Intercept                                     
    student_t(3, 0, 13.4)     sigma                                 0   
       source
      default
      default
 (vectorized)
 (vectorized)
      default
      default

mitchell.brm4 <- brm(mitchell.form,
  data = mitchell_sub,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 1000,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

mitchell.brm4 |>
  conditional_effects("Day") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm4 |>
  conditional_effects("Day") |>
  plot()

mitchell.brm4 |>
  conditional_effects("TOD") |>
  plot()

mitchell.brm4 |>
  conditional_effects("TOD") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

# mitchell.brm4 |> SUYR_prior_and_posterior()

9.0.1 stan plots

The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

mitchell.brm4$fit |> stan_trace()

The chains appear well mixed and very similar

stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

mitchell.brm4$fit |> stan_ac()

There is no evidence of auto-correlation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

mitchell.brm4$fit |> stan_rhat()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

mitchell.brm4$fit |> stan_ess()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

mitchell.brm4$fit |> stan_dens(separate_chains = TRUE)

9.0.2 pp check

mitchell.brm4 |> pp_check(type = "dens_overlay", ndraws = 100)

9.0.3 DHARMa residuals

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

## mitchell.resids <- make_brms_dharma_res(mitchell.brm4, integerResponse = FALSE)
## wrap_elements(~testUniformity(mitchell.resids)) +
##                wrap_elements(~plotResiduals(mitchell.resids, form = factor(rep(1, nrow(mitchell_sub))))) +
##                wrap_elements(~plotResiduals(mitchell.resids, quantreg = TRUE)) +
##                wrap_elements(~testDispersion(mitchell.resids))

## mitchell.resids |> testTemporalAutocorrelation(time = mitchell_sub$Day)

# mitchell.resids |> testTemporalAutocorrelation(time=unique(mitchell_sub$Day))
autocor_check(mitchell_sub, mitchell.brm3, variable = "Day", n.sim = 250)

variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram

Using  as id variables

`geom_smooth()` using formula = 'y ~ x'

autocor_check(mitchell_sub, mitchell.brm4, variable = "Day", n.sim = 250)

variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram

Using  as id variables
`geom_smooth()` using formula = 'y ~ x'

residuals(mitchell.brm3)[, "Estimate"] |> acf()

residuals(mitchell.brm4)[, "Estimate"] |> acf()

Conclusions:

the simulated residuals do suggest that there might be a dispersion issue
it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.

mitchell.brm3 |>
  augment() |>
  ## group_by(ID) |>
  reframe(ACF = as.numeric(acf(.resid, plot = FALSE)$acf)) |>
  ## group_by(ID) |>
  mutate(N = 1:n()) |>
  ## slice(1:12) |>
  ungroup() |>
  group_by(N) |>
  summarise(ACF = mean(ACF)) |>
  ggplot(aes(y = ACF, x = N)) +
  geom_hline(yintercept = 0) +
  geom_segment(aes(yend = 0, xend = N))

mitchell.brm4 |>
  augment() |>
  ## group_by(ID) |>
  reframe(ACF = as.numeric(acf(.resid, plot = FALSE)$acf)) |>
  ## group_by(ID) |>
  mutate(N = 1:n()) |>
  ## slice(1:12) |>
  ungroup() |>
  group_by(N) |>
  summarise(ACF = mean(ACF)) |>
  ggplot(aes(y = ACF, x = N)) +
  geom_hline(yintercept = 0) +
  geom_segment(aes(yend = 0, xend = N))

References

Mitchell, David J, Antoine M Dujon, Christa Beckmann, and Peter A Biro. 2019. “Temporal autocorrelation: a neglected factor in the study of behavioral repeatability and plasticity.” Behavioral Ecology 31 (1): 222–31. https://doi.org/10.1093/beheco/arz180.